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Abstract. Stars are gravitationally stabilized fusion reactors changing their chemical composition while trans- 
forming light atomic nuclei into heavy ones. The atomic nuclei are supposed to be in thermal equilibrium 
with the ambient plasma. The majority of reactions among nuclei leading to a nuclear transformation are 
inhibited by the necessity for the charged participants to tunnel through their mutual Coulomb barrier. As 
theoretical knowledge and experimental verification of nuclear cross sections increases it becomes possible to 
refine analytic representations for nuclear reaction rates. Over the years various approaches have been made 
to derive closed-form representations of thermonuclear reaction rates (Critchfield 1972, Haubold and John 
1978, Haubold, Mathai and Anderson 1987). They show that the reaction rate contains the astrophysical 
cross section factor and its derivatives which has to be determined experimentally, and an integral part of the 
thermonuclear reaction rate independent from experimental results which can be treated by closed-form rep- 
resentation techniques in terms of generalized hypergeometric functions. In this paper mathematical/statisti 
cal techniques for deriving closed-form representations of thermonuclear functions, particularly the four 
integrals 
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will be summarized and numerical results for them will be given. The separation of thermonuclear functions 
from thermonuclear reaction rates is our preferred result. The purpose of the paper is also to compare 
numerical results for approximate and closed-form representations of thermonuclear functions. This paper 
completes the work of Haubold, Mathai, and Anderson (1987). 
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1 Barrier penetration at astrophysical energies 

The majority of nuclear reactions of astrophysical interest are inhibited by the necessity for the charged 
participants to tunnel through their mutual Coulomb barrier. Nuclear processes such as a-decay and decay 
by emission of heavier nuclei are also mediated by penetration through a static, one-dimensional Coulomb 
potential barrier. Barrier penentration factors in nuclear reaction rates take into account the exponential 
nature of the tail of the nuclear potential. A great impact on the nature of the potential has also the 
inclusion of the electron screening of the reacting particles, which leads to potentials of the Yukawa type, 
which exhibits a change of the height and width of the barrier compared to the Coulomb type of potential 
(Fowler 1984). 

In order to extrapolate measured nuclear cross sections a(E) down to astrophysical energies, the nuclear 
cross section factor S(E) is introduced by 

a(E) = ^-exp{-2n V }, (1.1) 

where i] is the Sommerfeld parameter 

Z x Z 2 e 2 (nWZ x Ztf? 
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with Z the atomic charge and v [E] the asymptotic relative velocity [kinetic energy] of the reacting nuclei 
(Fowler 1984). Thus, the cross section is given by the product of the cross section factor to be determined 
experimentally, the square of the de Broglie wavelength due to quantum mechanics (~ -E -1 ), and the barrier 
penetration factor. The quantity exp{— 27T7/} takes exclusively s-wave transmission into account, describing 
penetration to the origin through a pure Coulomb potential. Nuclear reactions rates are extremely sensitive 
to the precise numerical value in the argument of this exponential factor. The inclusion of uncertainties in 
the shape of the nuclear potential and contributions from non s-wave transmission, respectively, are very 
important for deriving specific nuclear reaction rates but do not change the overall energy dependence of 
the nuclear cross-section given in (1.1). Actually, uncertainties in the shape of the nuclear potential tail and 
contributions from non s-wave terms are only important for heavy-ion reactions. In the following we are 
focusing on reaction rates of the proton capture type, i.e. small value of the reduced mass fi and small value 
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of the atomic charge product i?ii? 2 e 2 . The main uncertainty in (1.1) lies in the variation of the cross section 
factor S(E) with energy, which depends primarily on the value chosen for the radius at which formation of 
a compound nucleus between two interacting nuclei or nucleons occurs (Brown and Jarmie 1990). 

The separation of the barrier penetration factor in (1.1) is based on the solution of the Schrodinger 
equation for the Coulomb wave functions. Therefore the cross section <r(E) in (1.1) can be parametrized 
even more precisely by either expanding S(E) into a Taylor series about zero energy because of its slow 
energy dependence, 



S(E) = 5(0) 
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(1.3) 



where S(0) is the value of S(E) at zero energy, and S'(0) and S"(0) are the first and second derivatives of 
S(E) with respect to energy evaluated at E=0, respectively, or to elaborate on the action integral 7(ri,r 2 ) 
to include effects due to the shape of the nuclear potential and non s-wave contributions, where n and r 2 
are the inner turning point and outer turning point, respectively, where the reacting particles tunnel from 
n to r 2 in a Coulomb plus nuclear field (Smith, Kawano, and Malaney 1993). Then the barrier penetration 
factor in (1.1) can be expressed in terms of this action integral as 

T = exp{-2I(r 1 ,r 2 )}, (1.4) 

which simplifies for a Coulomb field and for n = to be / c (0,r 2 ) = ttij, where / c (0,r 2 ) is the sharp-cutoff 
Coulomb integral. If one takes into account non s-wave terms and does not confine to the sharp-cutoff 
approximation of the Coulomb integral in (1.4), the overall energy dependence of the nuclear cross section 
a(E) can be approximated by 

a(E) = + C 2 E^ + C 3 (C 4 + EW) + ...}, (1.5) 

where the leading term containing C\E~ X / 2 corresponds to the exponential term in (1.1); C 2 , C3 and C 4 arc 
energy independent nuclear constants (Rowley and Merchant 1991). 

Electron screening of reacting nuclei brings about a considerable enhancement of nuclear reactions, par- 
ticularly in high-Z matter. The Coulomb potential is modified by the presence of a polarising cloud of 
electrons surrounding the positive ions. The potential seen by a reacting nucleus is found to be narrower 
than the Coulomb potential and quantum-mechanical tunneling through the barrier becomes easier. The 



barrier penetration factor in (1.1) taking into account a screened potential can be written in terms of a 
screening parameter t, 

a { E) = —exp\-2*{-) ft(£; + f)1/a |, (1-6) 

where t = Z\Z2e 2 K and K denotes the Debye-Hiickel length. Screened nuclear reaction rates are extremely 
sensitive to the precise numerical value of the argument of the exponential factor in (1.6). 



2 Evolution towards the Maxwellian equilibrium distribution 

It is a major assumption in deriving nuclear reaction rates that the reacting nuclei are supposed to be in 
thermal equilibrium with the ambient plasma. This assumption can be justified by comparing the char- 
acteristic time for significant energy exchanges by Coulomb collisions with the characteristic time it takes 
the nuclear reaction to produce the final nucleus. Generally the Coulomb collision time is many orders of 
magnitude smaller than the time to produce the final nucleus which is the natural condition that the nuclei 
are in thermal equilibrium with the ambient plasma. Thus the velocity distribution function of nuclei is 
Maxwell-Boltzmannian. The state of the plasma at time t is described by the distribution function nf(v, t), 
where n is the constant particle number density, v is the velocity variable, and v =\ v |. Conservation of 
mass and energy imply that 

J d 3 vf(v,t) = l, 

dW/(M) = -, (2.1) 
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where T is the constant kinetic temperature and \i denotes the mass. In a gravitationally stabilized stellar 
fusion reactor, as t — > oo, f(v,t) tends to the Maxwell-Boltzmann distribution function, 

In a thermonuclear plasma, the reaction rate arises from an integral of the nuclear cross section (equations 
(1.1) or (1.6)), times velocity, times the Maxwell-Boltzmann distribution of velocities (2.2), 
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It is evident from (2.3) that the kernel of the integral consists of a product of the steeply falling Maxwcll- 
Boltzmann distribution (2.2) and the rapidly rising cross section (1.1) or (1.6) to produce a not quite 
symmetrical peak, commonly called the Gamow peak. This peak justifies the fact, that reaction rates are 
extremely sensitive to the precise numerical values in the arguments of the exponential factors exhibiting the 
exponential nature of the tail of the nuclear potential and the exponential nature of the tail of the velocity 
distribution function. 

The Maxwell-Boltzmann distribution is a solution of the general nonlinear Boltzmann equation which 
itself reveals as notoriously complicated. The system of particles here is considered to be an infinite, spatially 
homogeneous and isotropic gas containing a variety of nuclei. It is also assumed that only binary reactions 
need to be taken into account, so that the Boltzmann equation applies. Additionally the assumption is 
made that the nuclear reactions are isotropic, i.e., the cross section a is independent of the collision angle. 
Maxwell established that the low-order moments of the distribution function effectively relax toward their 
equilibrium values in just a few mean collision times. This corresponds, as discussed before, to the property 
that the low-energy part of the distribution attains Maxwell- Boltzmannian form in such a time interval. 
Nonlinear relaxation has been discussed by Kac (1955). 

On several occasions the question has been raised whether there may exist intermediate distributions that 
will evolve in such a way that the high- velocity tail of the respective velocity distribution will, at certain typ- 
ically high velocities and for certain time-intervals, display significant enhancement or depletion with respect 
to the steady-state Maxwell-Boltzmann distribution. Such a modification of the tail away from the Maxwell- 
Boltzmann distribution would significantly change the Gamow peak in (2.3) and subsequently would alter the 
respective reaction rates among nuclei in the plasma of the gravitationally stabilized stellar fusion reactor. A 
certain type of nonequilibrium distribution functions have been studied by Krook and Wu (1976, 1977), Tjon 
and Wu (1979), and Barnsley and Cornille (1981) by investigating solutions of the Boltzmann equation which 
approach an equilibrium distribution when t — > oo in a nonuniform fashion. This nonuniformity is due to 
the high velocity tail of the distribution and indicates that linearization techniques can not be fully justified 
for high velocities even when the state of the physical system is close to Maxwcll-Boltzmannian behavior. 
Their model considerations, while studying the relaxation of solutions of the Boltzmann equation towards 
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the steady-state Maxwcll-Boltzmann distribution, encourage the investigation of reaction rates containing 
a modified Maxwell-Boltzmann distribution. Having discussed the energy dependence of the nuclear cross 
section in Section 1 and Maxwell-Boltzmann distribution in Section 2, respectively, the following four inte- 
grals can be derived, representing thermonuclear functions for four quite different physical conditions. The 
standard case of the thermonuclear function contains the nuclear cross section (1.1), the energy dependent 
term of the Taylor series in (1.3), and the steady-state Maxwell-Boltzmann distribution function (2.2), 

r°° i 

h{z,v) d U / y v e-y e - z y~'dy (2.4) 
Jo 

where y = E/kT and z — 2n(fx/2kT) 1 ^ 2 ZiZ2e 2 /%. Considering dissipative collision processes in the ther- 
monuclear plasma cut off of the high energy tail of the Maxwell-Boltzmann distribution may occur, thus we 
write for (2.4), 

h(z,d,v) d = f [ y"e-ye-*y~ h dy, (2.5) 
Jo 

where d denotes a certain typically high energy. 

Accomodating screening effects in the standard thermonuclear function we have to use the nuclear cross 

section (1.6) and the steady-state Maxwcll-Boltzmann distribution function (2.2) which leads to 

r°° i 
I 3 (z,t,v)= f y»e-ye-^y +t y- 2 dy, (2.6) 
Jo 

where t is the electron screening parameter. 

Finally, if due to plasma effects a depletion of the Maxwell- Boltzmann distribution has to be taken into 
account, the thermonuclear function can be written in the follwing form 

h(z,6,b,v) = f / y v e-ye- b y S e- z y~*dy, (2.7) 
Jo 

where the parameter S exhibits the enhancement or reduction of the high-energy tail of the Maxwcll- 
Boltzmann distribution. 

In the following Sections mathematical/statistical techniques for deriving closed- form representations of 
the four thermonuclear functions (2.4) - (2.6) will be summarized, their asymptotic forms will be given and 
numerical results for both of them derived. 
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3 Mathematical preliminaries 

First of all, we need to recall the gamma function, defined for complex z by 

/•OC 

r(z) = / f-^dt^iz) > o. 

Jo 

The definition of T(z) can be extended to the entire complex plane where it is analytic except for simple 
poles at and the negative real integers. An important property we shall need is the multiplication formula 

T{mz) = (2n) 1 ^m mz -iT(z)T(z + 1) • • • T(z + ^), (3.1) 

which is valid for all z and all integers m > 1. 
Definition. The function 



" ' " "' x i f ii"' + s)n« =1 r(i - aj - s) 



- Z - s ds,z^0, (3.2) 



bl ,... bq J 2^ J L q =m+1 r(i - bj - a)n? =n+1 r( 0j + s) ' 

is called the G-function and is originally due to Meijer (cp. Mathai and Saxena 1973). Here, i = V^T; to, n, 
p, and q are integers with < n < p and < m < q. In (3.2), and throughout this paper, an empty product 
is interpreted as unity (similarly an empty sum as zero). The a^j = 1, . . . ,p and bj, j = 1, . . . ,q are complex 
numbers such that no pole of T(bj + s),j = 1, . . . , to coincides with any pole of T(l — a,j — s), j = 1, . . . , n. L 
is a contour separating the poles of T(bj + s),j = 1, . . . , m from the poles of T(l — aj — s), j = 1, . . . , n. At 
this point, it is not clear that the integral in (3.2) even exists. Conditions on the contour L and the various 
parameters must be imposed in order that the integral converges. These conditions, as well as properties of 
the G-function may be found in Luke (1969), chapter 5. However, for the G-functions encountered in this 
paper, it suffices to know that the integral in (3.2) is well-defined for all z ^ if 

(i) L is a loop beginning and ending at — oo and encircling all poles of T(bj + s), j = 1, . . . , m, once in the 
positive direction, but none of the poles of T(l — aj — s), j = 1, . . . , n, and 

(ii) q > 1 and p < q. 

Moreover, under these conditions the integral can be evaluated as a sum of residues at the poles of T(bj + 
s),j = 1, . . . ,TO. 
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One property that we will certainly require in the sequel is the asymptotic behaviour of G^(z) as 
\z\ — ► oo. From Luke (1969) page 179, we have 



G q '° 



( ^^-e'"^ z 6 as |z| -» oo, | arg^| < (a + e)ir - J, (3.3) 



where 
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and (70 = (1 - cr) + ^ - ^ «i 
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Definition. Let /(t) be a function defined for t > 0. Then 
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Jo 

/(i) = — / t-'M f {s)ds, (3.5) 
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is called a Mellin transform pair. (3.4) is called the Mellin transform, and (3.5) is the inversion formula. The 
transform normally exists only in the strip a < Sft(s) < (3, and the inversion contour must lie in this strip. 

Lemma 3.1 Let fi(t) and f2(t) be two functions with Mellin transforms Mf 1 (s) and Mf 2 (s). Then 

J o v- 1 f 1 {v)f2{-)dv^—J^M h (s)M f Mu- s ds. (3.6) 

Proof. Our proof is statistical. We suppose that f\{t) > 0, /2(i) > 0, J °° fi(t)dt < oo, and / °° f2(t)dt < oo 
(the application below will satisfy these criteria). By scaling if necessary, we can assume that fi(t) and 
f2(t) are density functions. Let X and Y be independent random variables having density functions fi(t) 
and f2(t) respectively. Then the left-hand side of (3.6) is the density function g(u) of the random variable 
U = XY. Let us look at the right-hand side. We have M fl (s) = E{X S ~ 1 ) and M h (s) = E(Y S ^ 1 ), and 
therefore 

M g (s) = EiU 3 - 1 ) = EiX-^Y 3 - 1 ) = E{X S ~ 1 )E{Y S - 1 ) = M h (s)M h {s). 

It follows that the right-hand side of (3.6) is ^ J L M g (s)u~ s ds 7 which is the formula for the inverse Mellin 
transform of M g (s). Thus the right-hand side is g{u) as well. 
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4 Representation of the four integrals in terms of G-functions 

Theorem 4.1 (Saxena (1960), Mathai and Haubold (1988)) For z > 0,p > 0, p < 0, and integers 
m, n > 1, we have 



p I t- np e- pt e~ zt ~ m dt = p np (2ir 
la 



n PC)-Tr\l( 2 - n - m ) 



m 2 n 2 
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(4.1) 
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Proof. Define f x {t) = i 1 "™^"* and f 2 (t) = e"* m for t > 0. Then the Mellin transforms are 

M fl (s)= t s - 1 f 1 (t)dt= t 1 - np+s - 1 e- t dt = T(l-np + s),n(l-np + s)>0, 
Jo Jo 



and 

M /2 ( S )= / 1f~ 1 f 2 {t)dt = / ^ 

Jo Jo 
Then by setting v = pt and u — z^p, and using the lemma, we have 
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e~ tm dt= -T(-s),$t(s) >0. 
n n 
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t -np e -pt e -zt m dt = p np J v -np 



• v- np e- v e-^ m dv=p np V 1 f 1 (v)f 2 (-)dv 
Jo Jo v 



= V — f M fl {s)M f2 {s)u- s ds = ^- [ T{l-np + s)-T(-s)u- s ds 
2m J L /lV ' i2K ' 2mJ L y F ' n K n ' 



rap 



T(l -np + ns')T(ms'){z m p n )- s ds 



2m J L , 

where we made a change of variable s = ns' . The G-fimction appearing on the right-hand side of (4.1) is 



(4.2) 



^m+n,0 I z"*p" 

\ o,i,...,^ ; i^i£,...,i!^ri£y 

= mk ns)n± + -)■■■ r(*=i + -)r(i=p* + «)... r(^ + ,) ~° ds. 

By the multiplication formula in (3.1), we have 

T(l - np + ns) = T(n[- -p + si) = (2^)^ n n ^- p+s ^-^T^—^- + s) ■ ■ ■ r(^— ^ + s). 

n n n 

Thus by applying the multiplication formula and (4.4) to (4.3), we get 



(4.3) 



(4.4) 
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j r(ms)T(l-np + ns){z m p n )- s ds. 



(4.5) 
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By comparing (4.2) and (4.5), we obtain (4.1). By setting m = 2, n = l,p = 1, and p = —u, we obtain 



Corollary 4.2 For z > and v > 0, we have 



Ii (*,!/)= / j/ I/ e-fe-^ 2 rf 2 / = 7r^G^ 



The proof of the following theorem is similar to that of theorem 4.1. 



0,1,1+"/ 



(4.6) 



Theorem 4.3 (Mathai and Haubold (1988)) For z > 0, d > 0, a > 0, and integers m, n > 1, we have 
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j-npg-atg-zt m dt = ^±f 2 ^^p. d l-np 
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p+l±l4-2_i, F l,...,„;2zl, 3 = l,..., m/ 



By setting m = 2, n = 1, a = 1, and p = — i/, we obtain 



Corollary 4.4 For z > 0, d > 0, and f > 0, we have 

h{z,d,v) = / 
Jo 



y» e -V e -*V 



The integral J3 may be worked out in terms of Ii and I2 as follows. We have 
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y+r+1,0 ± 



(4.7) 
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h(z,S,b,u) = f / ^e-^e-^ e-*y^dy= J/^V^J/ 
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In order to confidently exchange the summation and integral signs in (4.9), the quantity 
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(4.8) 



(4.9) 



11 



must be finite (by Fubini's theorem). Hence we expect the expansion in (4.9) may not be valid for large b 
and 6. (This was in fact borne out by later numerical computations.) 

To end this section, we use (3.3) to obtain asymptotic formulas for the four integrals. By a direct 
application of (3.3) to (4.6), (4.7), (4.8), and (4.9) (and some algebra in the case of the last three), wc obtain 



h(z,v) ~ 2 (J) 

h(z,d,u) ~ d v+1 e- d 

h(z,t,v) ~ 2(1)" 

h(z,6,b,v) ~ 2(|) f (^ 
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(4.10) 
(4.11) 
(4.12) 

(4.13) 



all as z — > oo. 



5 Series representations for the four integrals 

Series expressions for the four integrals can now be obtained by evaluating the G-functions using residue 
calculus. We will illustrate the method by doing this for the integral I\{z,u). This means that we have to 
evaluate the complex integral 



r»3>o ( z2 

G 0,3 t 



^//« r (5+ -Ml + " + «)(£) ^) 



o,|,i+^ 

As previously mentioned, the right-hand side will be the sum (i?i + i? 2 + -R3 below) of the residues of the 
integrand. We will assume that v is a non-negative integer (the analysis is slightly different otherwise, as 
seen in Mathai and Haubold (1988)). Then the poles of the gammas in the integrand of (5.1) are as follows: 

Poles of T(s): s = 0, -1, -2, . . . 

Poles ofr(i + .s): a = -|, -f, -§, . . . 

Poles of r(l + v + s); -v - 1, -v - 2, -v - 3, . . .. 
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Note that T(s) and T(l + v + s) have some poles in common. These will be poles of order two. Thus the 
poles 

s = 0, — 1, — 2, . . . , — i^are of order 1 each, 

s = -i, A A-- are of order leach, 
2 2 2 

s = — v — 1, — v — 2, — v — 3, . . . are of order 2 each. 



Using the facts that 



Hm„(,s + r)T(s) = T(a - r) = { .^[^ , r = 0, 1, 2, . . . ; T ( \ 

r! (1 — a) r \2 



7T 2 . 



where 

a(a + l)---(a + r- 1) ifr>l, 

1 if r = 0, 

we find that the sum of residues of the integrand at the poles s = 0, — 1, . . . , — u is 



(a) r 



r=0 



2\ - s 



E L T lr( 5- r)r(1+l '- r HT 



In exactly the same way, the sum of the residues at the poles s = 



1 _3 _5 
2' 2' 2' 



IS 



where 0F2 is the hypergeometric function defined by 

Finally, the sum of the residues at the poles s = — v — 1, — v — 2, . . . (each of order 2) is 



d 



r=0 



1 



(a + 1 + 1/ + ryT{s)T(- + s)T{l + v + s) — 



r=0 



log ( - ] + A r 
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where 



= [ l+i+... + i]- f 
+ 



± ~ 2 ~ ' r+v+1 



1/2 ~ 3/2 T 



37-21og2 



(l/2)+^+r_ 

(we take [1 + \ -\ h ±] to be zero if r = 0), 7 = 0.5772156649... is Euler's constant, and 



(5.2) 



_£>. r 



(_l)l+^+r r (_l _„) 



r!(r + 1/ + 1)!(|+ i/) r " 
By summing i? 1; i? 2 , and R3 and using (4.6) and (5.1), we obtain the following theorem. 

Theorem 5.1 Let v > be an integer. Then for z > 0, we have 

>2\5 / 3 ] 



(5.3) 



7lM =r(1 + ^S(IU^iv 4 



n . OO / 9 



r=0 



log 



A. 



where A r and -B r are given above in (5.2) and (5.3). The integral I 2 {z, d, v) can be treated in the same way 
with the following result. 

Theorem 5.2 Let z > 0, t > 0, and let v > be an integer. Then 



r=0 



+ £ 



-His) £ 

(£)' 



+ 



-/!(3) /(l/ + r +l_/) 
/ 2 \ "+r+l 

^ 



+r 



log 



4rf 



A 



where 



A = -27- 2 log 2- 



+ 



+ 



11 1 

I72 + 372 + " ' + (1/2) + v + r 



v + r + 1 

Since the integrals I3 and Z4 have been expressed in (4.8) and (4.9) in terms of I± and I 2 , similar 
expansions can and have been derived for ^3 and I4. However, the exact details will not be given here. 



6 Computations and Conclusion 

Numerical computations for the series expansions obtained above of the four integrals I\ , I 2 , I3, and J4 were 
made and compared to the corresponding approximations for large z in (4.10)-(4.13). The programming was 
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carried out in Pascal on a Macintosh II computer with a numerical coprocessor, for a wide range of parameter 
values. Some of the results are shown in figures 1 to 4. Obviously, the series computations fail for large 
values of z. Since much effort was made in optimizing the program for accuracy and countering problems 
of underflow and overflow, it is thought that this failure is a result of machine and compiler numerical 
accuracy. It is evident, however, that the missing portions of the "exact" curves can be replaced by the 
"approximate" curves. 

For the sake of comparison, the four integrals were computed for the same parameter values using the 
numerical integration routines in Mathematica (Wolfram 1991). The results were identical to the results of 
the previous paragraph, except that computations for larger values of z were possible. The results, together 
with corresponding approximations, are plotted in figures 5 to 8. 
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